{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://matplotlib.org/users/customizing.html\n",
    "# print(plt.style.available) # uncomment to print all styles\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=2)\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "mpl.rcParams['figure.figsize'] = (10.0, 8.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1) Density of Floating Point Numbers\n",
    "\n",
    "The following code snippet generates all possible floating point numbers in a floating point system and shows them in a plot to illustrate their density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "significand_bits = 2\n",
    "exponent_min = -4\n",
    "exponent_max = 4\n",
    "\n",
    "fp_numbers = []\n",
    "for exp in range(exponent_min, exponent_max+1):\n",
    "    for sbits in range(0, 2**significand_bits):\n",
    "        significand = 1 + sbits/2**significand_bits \n",
    "        fp_numbers.append(significand * 2**exp)\n",
    "        \n",
    "fp_numbers = np.array(fp_numbers)\n",
    "fp_numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,1))\n",
    "plt.plot(fp_numbers, np.ones_like(fp_numbers), \"o\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,1))\n",
    "plt.plot(fp_numbers, np.ones_like(fp_numbers), \"o\");\n",
    "plt.axis([0,1.6, 0,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2) Picking apart a floating point number (IEEE double precision)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Never mind the details of this function...\n",
    "\n",
    "def pretty_print_fp(x):\n",
    "    print(\"---------------------------------------------\")\n",
    "    print(\"Floating point structure for %r\" % x)\n",
    "    print(\"---------------------------------------------\")\n",
    "    import struct\n",
    "    s = struct.pack(\"d\", x)\n",
    "\n",
    "    def get_bit(i):\n",
    "        byte_nr, bit_nr = divmod(i, 8)\n",
    "        return int(bool(\n",
    "            s[byte_nr] & (1 << bit_nr)\n",
    "            ))\n",
    "\n",
    "    def get_bits(lsb, count):\n",
    "        return sum(get_bit(i+lsb)*2**i for i in range(count))\n",
    "\n",
    "    # https://en.wikipedia.org/wiki/Double_precision_floating-point_format\n",
    "\n",
    "    print(\"Sign bit (1:negative):\", get_bit(63))\n",
    "    \n",
    "    stored_exponent = get_bits(52, 11)\n",
    "    fraction = get_bits(0, 52)\n",
    "    \n",
    "    if stored_exponent == 0:\n",
    "        exponent = -1022\n",
    "        significand = fraction\n",
    "    else:\n",
    "        exponent = stored_exponent - 1023\n",
    "        significand = fraction + 2**52\n",
    " \n",
    "    print(\"Mathematical exponent (decimal): m = %d\" % exponent)\n",
    "    print(\"Shifted exponent value (decimal): c = %d\" % stored_exponent)\n",
    "    print(\"Stored exponent value (bin): c = {0:011b}\".format(stored_exponent) )\n",
    "\n",
    "    print(\"Fractional part of significand (bin): f = {0:052b}\".format(fraction))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pretty_print_fp(float(1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Things to try:\n",
    "\n",
    "* Twiddle the sign bit\n",
    "* 1,2,4,8\n",
    "* 0.5,0.25\n",
    "* $2^{\\pm 1023}$, $2^{\\pm 1024}$\n",
    "* `float(\"inf\")`, `float(\"-inf\")`\n",
    "* `float(\"nan\")`\n",
    "* $2^{\\pm 1040}$\n",
    "* `float(\"+0\")`, `float(\"-0\")`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin(1023)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pretty_print_fp(float(0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3) Checking the range of integers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discussed in class that we can obtain the range of integers that can be represented exactly as\n",
    "\n",
    "$$[1,2^{p}]$$\n",
    "\n",
    "So let's try this using Python (where double precision gives p=53)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This number is within the range\n",
    "2**52"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We add one to get the next integer \n",
    "2**52 + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All good so far, right? Let's check the integer:\n",
    "2**53"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If I add one, what do you expect us to get?\n",
    "type(2**53 + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Humm... not what you expected right? Why do you think this happened?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(2**53)\n",
    "print(2**53 + 1.0)\n",
    "print(2**53 + 2.0)\n",
    "print(2**53 + 3.0)\n",
    "print(2**53 + 4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integer type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to find the maximum integer (using 64-bit signed integer)\n",
    "# uses only 63 bits, since one is reserved for the sign\n",
    "maxint = 0\n",
    "for i in range(63):\n",
    "     maxint += 2**i\n",
    "maxint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4) What will happen when we run these code snippets?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A) it won't stop (infinite loop)\n",
    "\n",
    "B) it will stop when b reaches underflow\n",
    "\n",
    "C) it will stop when b reaches machine epsilon\n",
    "\n",
    "D) none of the above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 10**4\n",
    "b = 1.0\n",
    "i=0\n",
    "while (a + b > a):\n",
    "    #print(\"{0:d}, {1:1.16f}\". format(i, b)) \n",
    "    print(i,b)\n",
    "    b = b / 2.0\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A) it won't stop (infinite loop)\n",
    "\n",
    "B) it will stop when b reaches underflow\n",
    "\n",
    "C) it will stop when b reaches machine epsilon\n",
    "\n",
    "D) none of the above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1.0\n",
    "while a > 0.0: \n",
    "    a = a / 2.0\n",
    "    print(a)\n",
    "    #print(\"% .16e\"% (a)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A) it won't stop (infinite loop)\n",
    "\n",
    "B) it will stop when b reaches underflow\n",
    "\n",
    "C) it will stop when b reaches machine epsilon\n",
    "\n",
    "D) none of the above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1.0\n",
    "while a != np.inf:\n",
    "    a = a * 2.0\n",
    "    print(a)\n",
    "    #print(\"% .16e\"% (a)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's write the following number (not exacly the definition of OFL, but close)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "float(2**1023)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try the UFL definition\n",
    "2**(1023+1)*(1-2**(-53))\n",
    "\n",
    "# Can we make this work?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
